from pyexpat import model
import numpy as np
import pandas as pd

# TSA from Statsmodels
import statsmodels.api as sm
import statsmodels.formula.api as smf
import statsmodels.tsa.api as smt

# Display and Plotting
import matplotlib.pylab as plt
import seaborn as sns


def arima1():
    import pandas 
    
    # 读取数据，指定日期为索引列
    
    data = pandas.read_csv(
        'D:\\DATA\\pycase\\number2\\9.3\\Data.csv' ,
        index_col='日期'
    )
    
    # 绘图过程中
    
    import  matplotlib.pyplot as plt
    
    # 用来正常显示中文标签
    
    plt.rcParams['font.sans-serif']=['SimHei']
    
    # 用来正常显示负号
    
    plt.rcParams['axes.unicode_minus'] = False 
    
    # 查看趋势图
    data.plot() #有增长趋势，不平稳
    
    
    # 附加：查看自相关系数合片自相关系数（查分之后），可以用于平稳性的检测，也可用于定阶系数预估
    
    #自相关图（）
    
    from statsmodels.graphics.tsaplots import plot_acf
    
    plot_acf(data).show() #自相关图既不是拖尾也不是截尾。以上的图的自相关是一个三角对称的形式，这种趋势是单调趋势的典型图形，说明这个序列不是平稳序列
    
    

    # 1 平稳性检测
    
    from statsmodels.tsa.stattools import adfuller as ADF
    
    
    def tagADF(t):
        result = pandas.DataFrame(index=[
                "Test Statistic Value", "p-value", "Lags Used", 
                "Number of Observations Used", 
                "Critical Value(1%)", "Critical Value(5%)", "Critical Value(10%)"
            ], columns=['销量']
        )
        result['销量']['Test Statistic Value'] = t[0]
        result['销量']['p-value'] = t[1]
        result['销量']['Lags Used'] = t[2]
        result['销量']['Number of Observations Used'] = t[3]
        result['销量']['Critical Value(1%)'] = t[4]['1%']
        result['销量']['Critical Value(5%)'] = t[4]['5%']
        result['销量']['Critical Value(10%)'] = t[4]['10%']
        return result
    
    
    print('原始序列的ADF检验结果为:',tagADF(ADF(data[u'销量'])))  # 添加标签后展现
    
    # 平稳判断：得到统计量大于三个置信度(1%,5%,10%)临界统计值，p值显著大于0.05，该序列为非平稳序列。
    # 备注：得到的统计量显著小于3个置信度（1%，5%，10%）的临界统计值时，为平稳 此时p值接近于0 此处不为0，尝试增加数据量，原数据太少
    
    # 2 进行数据差分，一般一阶差分就可以
    
    D_data = data.diff(1).dropna()
    D_data.columns = [u'销量差分']
    
    #差分图趋势查看
    
    D_data.plot() 
    plt.show()
    
    # 附加：查看自相关系数合片自相关系数（查分之后），可以用于平稳性的检测，也可用于定阶系数预估
    
    #自相关图
    
    plot_acf(D_data).show()
    
    plt.show()
    
    #偏自相关图
    
    from statsmodels.graphics.tsaplots import plot_pacf
    
    plot_pacf(D_data).show()
    
    # 3 平稳性检测
    
    print(u'差分序列的ADF检验结果为：', tagADF(ADF(D_data[u'销量差分']))) 
    
    # 解释：Test Statistic Value值小于两个水平值，p值显著小于0.05，一阶差分后序列为平稳序列。
    
    # 4 白噪声检验
    from statsmodels.stats.diagnostic import acorr_ljungbox
    
    #返回统计量和p值
    
    print(u'差分序列的白噪声检验结果为：', acorr_ljungbox(D_data, lags=1))  # 分别为stat值（统计量）和P值
    
    # P值小于0.05，所以一阶差分后的序列为平稳非白噪声序列。
    
    
    # 5 p，q定阶
    
    from statsmodels.tsa.arima_model import ARIMA
    
    #一般阶数不超过length/10
    
    pmax = int(len(D_data)/10) 
    
    
    #一般阶数不超过length/10
    
    qmax = int(len(D_data)/10) 
    
    #bic矩阵
    
    bic_matrix = [] 
    for p in range(pmax+1):
        tmp = []
    for q in range(qmax+1):
    #存在部分报错，所以用try来跳过报错。
        try: 
            tmp.append(ARIMA(data, (p,1,q)).fit().bic)
        except:
            tmp.append(None)
    bic_matrix.append(tmp)
    
    #从中可以找出最小值
    
    bic_matrix = pandas.DataFrame(bic_matrix) 
    
    #先用stack展平，然后用idxmin找出最小值位置。
    
    p,q = bic_matrix.stack().idxmin() 
    
    
    
    print(u'BIC最小的p值和q值为:%s、%s' %(p,q))
    # 取BIC信息量达到最小的模型阶数，结果p为0，q为1，定阶完成。
    
    # 6 建立模型和预测
    
    model = ARIMA(data, (p,1,q)).fit() 
    
    #给出一份模型报告
    
    model.summary2() 
    
    #作为期5天的预测，返回预测结果、标准误差、置信区间。
    
    model.forecast(5)
 

def arima2():
    ts_df = pd.read_excel(r"./data/arima.xlsx")
    n_sample = ts_df.shape[0]

    n_train = int(0.95*n_sample) + 1
    n_forecast = n_sample - n_train
    ts_train = ts_df.iloc[:n_train]['value']
    ts_test = ts_df.iloc[n_train:]['value']

    # 用来画四格图的方法
    # ts_ax, acf_ax, pacf_ax, plt = tsplot(ts_train, title='A Given Training Series', lags=20)
    # plt.show()

    # 用来查看参数的选取
    # get_best_parameters_arima2(train=ts_train)

    # Model Estimation
    # Fit the model

    arima200 = sm.tsa.SARIMAX(ts_train, order=(2,0,0))  # 这个元组就是(p,d,q)
    model_results = arima200.fit() 
    #残差分析 正态分布 QQ图线性
    model_results.plot_diagnostics(figsize=(16, 12));#statsmodels库
    plt.show()
    
def get_best_parameters_arima2(train):
    import itertools

    ts_train = train
    #当多组值都不符合时，遍历多组值，得出最好的值
    p_min = 0
    d_min = 0
    q_min = 0
    p_max = 4
    d_max = 0
    q_max = 4
    
    # Initialize a DataFrame to store the results
    results_bic = pd.DataFrame(index=['AR{}'.format(i) for i in range(p_min,p_max+1)],
                            columns=['MA{}'.format(i) for i in range(q_min,q_max+1)])
    
    for p,d,q in itertools.product(range(p_min,p_max+1),
                                range(d_min,d_max+1),
                                range(q_min,q_max+1)):
        if p==0 and d==0 and q==0:
            results_bic.loc['AR{}'.format(p), 'MA{}'.format(q)] = np.nan
            continue
        
        try:
            model = sm.tsa.SARIMAX(ts_train, order=(p, d, q),
                                #enforce_stationarity=False,
                                #enforce_invertibility=False,
                                )
            results = model.fit()
            results_bic.loc['AR{}'.format(p), 'MA{}'.format(q)] = results.bic
        except:
            continue
    results_bic = results_bic[results_bic.columns].astype(float)

    # 画图
    fig, ax = plt.subplots(figsize=(10, 8))
    ax = sns.heatmap(results_bic,
                    mask=results_bic.isnull(),
                    ax=ax,
                    annot=True,
                    fmt='.2f',
                    )
    ax.set_title('BIC')
    plt.show()


def tsplot(y, lags=None, title='', figsize=(14, 8)):
    
    fig = plt.figure(figsize=figsize)
    layout = (2, 2)
    ts_ax   = plt.subplot2grid(layout, (0, 0))
    hist_ax = plt.subplot2grid(layout, (0, 1))
    acf_ax  = plt.subplot2grid(layout, (1, 0))
    pacf_ax = plt.subplot2grid(layout, (1, 1))
    y.plot(ax=ts_ax)
    ts_ax.set_title(title)
    y.plot(ax=hist_ax, kind='hist', bins=25)
    hist_ax.set_title('Histogram')
    smt.graphics.plot_acf(y, lags=lags, ax=acf_ax)
    smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax)
    [ax.set_xlim(0) for ax in [acf_ax, pacf_ax]]
    sns.despine()
    plt.tight_layout()
    return ts_ax, acf_ax, pacf_ax, plt
    


if __name__ == '__main__':
    arima2()